La transcriptomique est l’étude des ARNm qui sont transcrits à partir d’un génome. Elle permet de quantifier le niveau d’expression d’un gène transcrit en ARNm qui est ensuite traduit en protéine. La mesure de l’expression des gènes peut se faire par la méthode de la PCR quantitative (quantitative Polymerase Chain Reaction ou qPCR). On obtient ainsi une mesure relative de l’expression d’un gène par rapport à un contrôle interne (normalisation par rapport à un/plusieurs gène(s) de référence).
Les données à analyser correspondent à des niveaux d’expression de micro-ARN (miR), soit des courtes chaines d’ARN (en moyenne 22 nucléotides) qui sont des régulateurs traductionnels: en s’appariant à un ARNm complémentaire, il va le dégrader ou réprimer sa traduction en protéine.
Des données d’expression de 19 miRs mesurées par qPCR sont disponibles et l’on vous demande de les analyser à l’aide de certaines méthodes multivariées vues durant les sessions théoriques.
Vous disposez pour cela du jeu de données data_qPCR.csv qui comporte les variables “Groupe” et “Identifiant” caractérisant chaque observation, ainsi que les données d’expression de chaque miR.
Avant de faire ces analyses, créez des matrices vides pour stocker efficacement les différentes sorties suivantes pour chaque classifieur (le classifieur est mis en ligne dans la matrice): * Une matrice reprenant le RMSE, RMSE-CV et la dimensionalité du modèle * Une matrice des coefficients de chaque méthode (sauf pour le Random Forest) * Une matrice des valeurs prédites par cross-validation pour chaque observation
Conseil: créez une fonction de cross-validation de type k-fold pour calculer plus facilement le RMSE-CV.
# chemin d'accès
data.path <- "/Users/bgovaerts/Dropbox/Partage_OMICS_COURSE/Data_stat2340/Data_qPCR"
# lecture du data frame
dataqPCR <- read.csv(file.path(data.path,"data_qPCR.csv"),sep=';',header=T)
cat("dimensions \n")## dimensions
dim(dataqPCR)## [1] 120 40
cat("colnames \n")## colnames
colnames(dataqPCR)## [1] "Groupe" "Identifiant" "miR_001_5p" "miR_002_5p"
## [5] "miR_003_3p" "miR_004_5p" "miR_005_5p" "miR_006_5p"
## [9] "miR_006_3p" "miR_007_5p" "miR_008_5p" "miR_009_5p"
## [13] "miR_010_3p" "miR_011_3p" "miR_012_5p" "miR_013_5p"
## [17] "let_01_5p" "let_02_5p" "miR_014_5p" "miR_015"
## [21] "miR_016" "Log_miR_001_5p" "log_miR_002_5p" "log_miR_003_3p"
## [25] "log_miR_004_5p" "log_miR_005_5p" "log_miR_006_5p" "log_miR_006_3p"
## [29] "log_miR_007_5p" "log_miR_008_5p" "log_miR_009_5p" "log_miR_010_3p"
## [33] "log_miR_011_3p" "log_miR_012_5p" "log_miR_013_5p" "log_let_01_5p"
## [37] "log_let_02_5p" "log_miR_014_5p" "log_miR_015" "log_miR_016"
# matice d'expression
X <- as.matrix(dataqPCR[,3:21]) # matrice d'expression
rownames(X) <- dataqPCR[,2]
obsnames=rownames(X)
X <- X[,sort(colnames(X))]
miR_names <- colnames(X)
Xcat=T # indique que les variables sont des catégories bien spécifiques et pas des longueurs d'onde par exemple. Cela a son importance dans les graphiques de coefficients et loadings
Xlog <- as.matrix(dataqPCR[,22:40])
rownames(Xlog) <- dataqPCR[,2]
Xlog <- Xlog[,sort(colnames(Xlog))]
n = dim(X)[1]
m = dim(X)[2]
classY <- rep(0, n)
classY[dataqPCR[,1]=="malade"] = 1
cat("class factor levels \n")## class factor levels
unique(classY)## [1] 1 0
table(classY, dataqPCR[,1])##
## classY malade temoin
## 0 0 57
## 1 63 0
# modification de matrices
# centrage miRs
XC <- scale(Xlog, center = T, scale = F)
XC <- data.frame(XC)
# On centre Y pour le PLS et l'O-PLS
YC <- scale(classY, center = T, scale = F)Nous allons ajuster 6 types de modèles (regression logistique SW, PLS, O-PLS et logisitique ridge et logistique LASSO, random forest) afin de prédire la classe en fonction des MIR.
Pour chacun des modèles, nous allons stocker le RMSE classique (training dataset = test dataset) ainsi que le RMSE obtenu par cross-validation k-fold. Pour certains modèles, nous allons également stocker les valeurs des coefficients.
RMSE <- function(y_true, y_fit) {
n <- length(y_true)
sqrt((1/n)*sum((y_true-y_fit)^2))
}
##################
# kfoldCV
##################
# fonction qui fait de la k-fold CV sur une fonction de régression prédéfinie prédisant y (FUN)
kfoldCV=function(X, Y, k, FUN = ypred_FUN, seed=100, ...) {
# Inputs
# X, la matrice des régresseurs
# Y, le vecteur des réponses
# k, le nombre de segments et
# FUN : la fonction de régression qui prédit Y
# Outputs :
# ypred.cv, les prédictions par CV
# 1. definition de k segments aleatoires
n=dim(X)[1]
set.seed <- seed
indx <- sample(1:n,n)
nseg <- k
sn <- ceiling(n/k) # taille des segments
# 2. boucle sur les sous-echantillons
# pour calculer y-fit par cv
yfit.cv <- numeric()
for (i in 1:nseg) {
if(i < nseg) {indseg <- indx[(i-1)*sn + (1 : sn)]}
else {indseg <- indx[((i-1) * sn + 1) : n]}
Y.train <- Y[-indseg]
X.train <- X[-indseg,]
X.valid <- X[indseg,]
yfitcv <- FUN(X = X.train, Y = Y.train, XV = X.valid, ...)
yfit.cv[indseg] <- yfitcv
}
# 3. Sortie de la fonction
return(yfit.cv)
}# Définition de paramètres de modélisation
# nombre de segments en cross validation
k <- 10
# Nombre de composantes max
NCmax=10
# Liste des modèles
modnames <- c('Logistique-SW','PLSR','O-PLS','Logistique-RIDGE','Logistique-LASSO','Random Forest')
nmodels <- length(modnames) # nombre de modèles évalués
# Preparation des matrices de résultats
# Matrice des RMSE
RMSEMat <- matrix(nrow = nmodels, ncol = 3)
colnames(RMSEMat) <- c('RMSE-train','RMSE-CV','Taille du modèle')
rownames(RMSEMat) <- modnames
# Matrice des coefficients (pour tous les modèles sauf RF)
CoefficientMat <- matrix(0, ncol = m, nrow = nmodels-1)
colnames(CoefficientMat) <- colnames(XC)[1:m]
rownames(CoefficientMat) <- modnames[-nmodels]
# Matrice des prédictions
PredMat <- matrix(nrow = n, ncol = nmodels)
colnames(PredMat) <- modnames
# Tout est mis dans une liste qui sera remis à jour modèle par modèle
listr=list(CoefficientMat=CoefficientMat,RMSEMat=RMSEMat,PredMat=PredMat)# création d'une fonction pour stocker et imprimer les résultats de chaque modèle
printresmod=function(i,YC,yfit,yfitCV,coef,parmod,listr,pcoef=T)
{
# i = N° du modèle
# YC : réponses observée
# yfit : réponses prédites
# yfitCV : réponses prédites par CV
# coef : vecteur des coefficients
# parmod paramètre de dimension ou de lissage du modèle
# pcoef : T/F pour indiquer si des coefficient sont dispo pour le modèle (pas le cas pour les RF)
par(mfrow=c(1,2))
plot(YC,yfit,pch=20,main="Y-Yfit")
plot(YC,yfitCV,pch=20,main="Y-YfitCV",ylab="yfitCV")
# calcul des RMSE
n=length(YC)
listr$RMSEMat[i,1]=sqrt(sum((YC-yfit)^2)/n)
listr$RMSEMat[i,2] =sqrt(sum((YC-yfitCV)^2)/n)
listr$RMSEMat[i,3] = parmod
# Récupération des coefficients et des réponses
if(pcoef==T) listr$CoefficientMat[i,]=coef
listr$PredMat[,i]=yfit
# affichage de résultats
pander(listr$RMSEMat[i,])
par(mfrow=c(1,1))
if(pcoef==T){
if(Xcat==F){plot(Xval,listr$CoefficientMat[i,],type="l",ylab="Coefficients",main=paste("Coefficients ",dimnames(coefficient)[[1]][i]),xlab=Xxlab)
abline(h=0)}
if(Xcat==T){dotchart(coef,labels=colnames(coef),
xlab="Coefficients",main=paste("Coefficients - ",modnames[i])) }
}
return(listr)
}# Ajustement du modèle de regression stepwise forward
null <- glm(classY ~ 1, data=XC,family="binomial") # modele de depart
full <- glm(classY ~ ., data=XC,family="binomial") # modele avec toutes les variables
fit.sw <- step(null, scope=list(lower=null, upper=full), direction="both",trace=0) # critère d'opimisation= AIC
pander(summary.glm(fit.sw))| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | -0.5942 | 0.3373 | -1.761 | 0.07817 |
| log_miR_002_5p | -12.62 | 3.313 | -3.81 | 0.0001387 |
| log_miR_009_5p | 13.01 | 3.304 | 3.939 | 8.197e-05 |
| log_miR_011_3p | 2.022 | 0.9884 | 2.046 | 0.04076 |
| log_miR_010_3p | -1.175 | 0.6565 | -1.79 | 0.07349 |
(Dispersion parameter for binomial family taken to be 1 )
| Null deviance: | 166.1 on 119 degrees of freedom |
| Residual deviance: | 104.1 on 115 degrees of freedom |
yfit.train <- fit.sw$fitted.values
namessw=names(coefficients(fit.sw))
coef <- CoefficientMat[1,]
coef[namessw[-1]]=coefficients(fit.sw)[-1]
nc <- length(coef)-1
ypred_FUN <- function(X,Y,XV,...){
null <- glm(Y ~ 1, data = X,family = "binomial")
full <- glm(Y ~ ., data = X,family = "binomial")
fit.sw <- step(null, scope = list(lower = null,
upper = full),
direction="both",trace=0)
coef <- coefficients(fit.sw)
XV <- cbind(1, XV[,names(coef)[-1]])
yfit.cv <- as.matrix(XV) %*% coef
yfit.cv <- exp(yfit.cv)/(1+exp(yfit.cv))
return(yfit.cv)
}
yfit.cv <- kfoldCV(X = XC, Y = classY, k = k, FUN = ypred_FUN)
# Sauvegardes et impression de résultats
listr=printresmod(1,classY,yfit.train,yfit.cv,coef,nc,listr)PLSR par la librairie pls. On estime directement les modèles par cross validation afin de trouver le bon nombre de composantes. Ensuite on reestime ce modèle là.
# Recherche du nombre de composantes optimales par CV
## définition de la fonction d'estimation du modèle
ypred_FUN <- function(X, Y, XV, ...){
fit.plsr <- plsr(Y ~ . , data = X, ...)
coef <- coefficients(fit.plsr)
yfit.train <- as.matrix(XV) %*% coef
return(yfit.train)
}
## Boucle sur le nombre de composantes possibles: 1 à NCmax
RMSE.cv <- numeric()
yfit.cv <- list()
for(i in 1:NCmax) {
res.cv <- kfoldCV(X = XC, Y = YC, k = k,
FUN = ypred_FUN,
ncomp = i)
yfit.cv[[i]] <- res.cv
RMSE.cv[i] <- RMSE(y_true = classY, y_fit = yfit.cv[[i]])
}
# Recherche du nombre de composantes à garder
plot(1:10, RMSE.cv, type = "o", main = "RMSE cross-validation")
ncomp.opt = which.min(RMSE.cv)
ncomp.opt = max(2, ncomp.opt)
abline(v = ncomp.opt, col = 2, lty = 2)
RMSE.cv <- RMSE.cv[ncomp.opt]
yfit.cv <- yfit.cv[[ncomp.opt]] + mean(classY)
# Ajustement du modèle avec le nombre optimal de composantes
fit.plsr <- plsr(YC ~ . , data = data.frame(XC),
ncomp = ncomp.opt)
print(summary(fit.plsr))## Data: X dimension: 120 19
## Y dimension: 120 1
## Fit method: kernelpls
## Number of components considered: 2
## TRAINING: % variance explained
## 1 comps 2 comps
## X 36.76 47.46
## YC 26.76 33.89
## NULL
yfit.train <- fit.plsr$fitted.values[, , ncomp.opt] + mean(classY)
coef <- coefficients(fit.plsr)[, , ]
nc <- ncomp.opt
# Sauvegardes et impression de résultats
listr=printresmod(2,classY,yfit.train,yfit.cv,coef,nc,listr)Représentation graphique des scores pour les 2 premières composantes avec la librairie MBX
par(mfrow=c(1,1))
# graphe des scores avec la librairie MBXUCL
ScatterPlot(fit.plsr$scores[,1],fit.plsr$scores[,2], createWindow=FALSE, points_labs =obsnames, main = paste("PLS score plot for PC1 and PC2"), color=classY, pch=classY,xlab="PC1",ylab="PC2")loadings <- fit.plsr$loadings
for (i in 1:2) {
plot(loadings[,i], type="h", xaxt="n", xlab="",
ylab = "Loading", main = paste0("Loading ", i), lwd = 3)
axis(side = 1, at = 1:length(miR_names), labels = miR_names, las=2, cex.axis=0.7)
}fit.opls <- opls(x = Xlog, y = classY, predI = 1,
orthoI = 4, scaleC = "center",
printL = TRUE, plotL = FALSE)
# Recherche du nombre de composantes optimales par CV
## définition de la fonction d'estimation du modèle
ypred_FUN=function(X,Y,XV,...) {
fit.opls <- opls(x = X, y = Y, predI = 1,
scaleC = "center" , ...)
ropls::predict(fit.opls, XV)
}
## Boucle sur le nombre de composantes possibles
RMSE.cv <- numeric()
yfit.cv <- list()
for(i in 1:NCmax) {
rescv <- kfoldCV(Xlog, classY, k = k, FUN = ypred_FUN,
orthoI = i, printL = FALSE,
plotL = FALSE)
yfit.cv[[i]] <- rescv
RMSE.cv[i] <- RMSE(y_true = classY, y_fit = yfit.cv[[i]])
}
# Recherche du nombre de composantes à garder
plot(RMSE.cv, type = "o", main = "RMSE cross-validation")
ncomp.opt <- which.min(RMSE.cv)
ncomp.opt <- max(2, ncomp.opt)
abline(v = ncomp.opt, col = 2, lty = 2)
yfit.cv <- yfit.cv[[ncomp.opt]]
# Ajustement du modèle avec le nombre optimal de composantes
fit.opls <- opls(x = Xlog, y = classY, predI = 1, orthoI = ncomp.opt,
scaleC = "center", printL = TRUE, plotL = FALSE)
yfit.train <- ropls::predict(fit.opls, X)
coef <- fit.opls@weightMN %*% fit.opls@cMN
nc <- ncomp.opt
# Sauvegardes et impression de résultats
listr=printresmod(3,classY,yfit.train,yfit.cv,coef,nc,listr)Représentation graphique des scores OPLS
scores_p <- fit.opls@scoreMN
scores_o <- fit.opls@orthoScoreMN
par(mfrow=c(1,1))
# graphe des scores avec la librairie MBXUCL
ScatterPlot(scores_p[,1], scores_o[,1], createWindow=FALSE, points_labs =obsnames, main = paste("O-PLS score plot for PCP and PCO1"), color=classY, pch=classY,xlab="PCPredictive",ylab="PCO1")loadings_p <- fit.opls@loadingMN
loadings_o <- fit.opls@orthoLoadingMN
plot(loadings_p, type="h", xaxt="n", xlab="", ylab = "Loading",
main ="Loading prédictif", lwd = 3)
axis(side = 1, at = 1:length(miR_names), labels = miR_names, las=2, cex.axis=0.7)
plot(loadings_o[,1], type="h", xaxt="n", xlab="", ylab = "Loading",
main ="Loading ortho 1", lwd = 3)
axis(side = 1, at = 1:length(miR_names), labels = miR_names, las=2, cex.axis=0.7)# Ajustement du modèle de regression pénalisé avec toutes les données
# lambda2 => RIDGE
fit.ridge <- penalized(response = classY,
penalized = XC,
model = "logistic",
lambda1 = 0,
lambda2 = 0.1,
trace = FALSE)
# Optimisation du paramètre lambda2 par cross-validation
lambda2 <- optL2(response = classY, penalized = XC,
lambda1 = 0, model = "logistic", trace=0)$lambda
# Regression logistique ridge avec lambda optimisé
fit.ridge <- penalized(response = classY,penalized = XC, model = "logistic",
lambda1 = 0, lambda2 = lambda2, trace = FALSE)
yfit.train = attr(fit.ridge,'fitted')
coef <- attr(fit.ridge,'penalized')
nc = sum(!(coef)==0)
# Calcul du yfit.cv
# fonction pour obtenir les y prédits
ypred_FUN = function(X, Y, XV, ...){
fit.penalized <- penalized(response = Y, penalized = X,
model = "logistic", trace=FALSE, ...)
penalizedcoeff <- attr(fit.penalized,'penalized')
intercept <- attr(fit.penalized,'unpenalized')
coef <- c(intercept, penalizedcoeff)
yfit.cv=as.matrix(cbind(1,XV))%*%coef
yfit.cv = exp(yfit.cv)/(1+exp(yfit.cv))
return(yfit.cv)
}
yfit.cv <- kfoldCV(X, classY, k = k, FUN = ypred_FUN, lambda1 = 0,
lambda2 = lambda2)
# Sauvegardes et impression de résultats
listr=printresmod(4,classY,yfit.train,yfit.cv,coef,nc,listr)# Ajustement du modèle de regression pénalisé avec toutes les données
# lambda1=> LASSO
fit.lasso <- penalized(response = classY, penalized = XC, model = "logistic",
lambda1 = 0.1, lambda2 = 0, trace = FALSE)
# Optimisation du paramètre lambda2 par cross-validation
lambda1 <- optL1(response = classY, penalized = XC,
lambda2 = 0, model = "logistic",
trace=0)$lambda
# Regression logistique ridge avec lambda optimisé
fit.lasso <- penalized(response = classY, penalized = XC, model = "logistic",
lambda2 = 0, lambda1 = lambda1, trace = FALSE)
yfit.train = attr(fit.lasso,'fitted')
coef <- attr(fit.lasso,'penalized')
nc = sum(!(coef)==0)
# Calcul du yfit.cv
# fonction pour obtenir les y prédits (c'est la m^ême que la précédente)
ypred_FUN = function(X, Y, XV, ...){
fit.penalized <- penalized(response = Y, penalized = X,
model = "logistic", trace=FALSE, ...)
penalizedcoeff <- attr(fit.penalized,'penalized')
intercept <- attr(fit.penalized,'unpenalized')
coef <- c(intercept, penalizedcoeff)
yfit.cv=as.matrix(cbind(1,XV))%*%coef
yfit.cv = exp(yfit.cv)/(1+exp(yfit.cv))
return(yfit.cv)
}
yfit.cv <- kfoldCV(X, classY, k = k, FUN = ypred_FUN,
lambda1 = lambda1, lambda2 = 0)
# Sauvegardes et impression de résultats
listr=printresmod(5,classY,yfit.train,yfit.cv,coef,nc,listr)fit.rf <- randomForest(formula=as.factor(classY) ~ ., data=Xlog, ntree=500, mtry=3)
# Recherche du nombre de composantes optimales par CV
## définition de la fonction d'estimation du modèle
ypred_FUN <- function(X,Y,XV,...) {
fit.rf <- randomForest(formula=as.factor(Y) ~ ., data = X, ntree = 500, ...)
yfit.cv <- predict(fit.rf, newdata = XV, type = 'class')
yfit.cv
}
## Boucle sur le nombre de composantes possibles
RMSE.cv <- numeric()
yfit.cv <- list()
j <- 1
mtry_tested <- seq(3,19, 2)
for(i in mtry_tested) {
rescv <- kfoldCV(Xlog, classY, k = k, FUN = ypred_FUN, mtry = i)
yfit.cv[[j]] <- rescv
yfit <- as.numeric(yfit.cv[[j]]) - 1
RMSE.cv[j] <- RMSE(y_true = classY, y_fit = yfit)
j <- j + 1
}
# Recherche du nombre de composantes à garder
plot(mtry_tested,RMSE.cv, type = "o", main = "RMSE cross-validation")
mtry.opt <- mtry_tested[which.min(RMSE.cv)]
abline(v = mtry.opt, col = 2, lty = 2)
yfit.cv <- as.numeric(yfit.cv[[which.min(RMSE.cv)]])-1
# Ajustement du modèle avec le nombre optimal de mtry
fit.rf <- randomForest(formula=as.factor(classY)
~ ., data = Xlog,
ntree=500, mtry=mtry.opt)
yfit.train <- predict(fit.rf, newdata = Xlog, type='class')
yfit.train <- as.numeric(yfit.train)-1
nc = m
# Sauvegardes et impression de résultats
listr=printresmod(6,classY,yfit.train,yfit.cv,coef,nc,listr,pcoef=F)cat(c("Taille des segments", k))## Taille des segments 10
pander(listr$RMSEMat)| RMSE-train | RMSE-CV | Taille du modèle | |
|---|---|---|---|
| Logistique-SW | 0.3829 | 0.443 | 18 |
| PLSR | 0.406 | 0.4301 | 2 |
| O-PLS | 0.5736 | 0.4473 | 2 |
| Logistique-RIDGE | 0.399 | 0.4442 | 19 |
| Logistique-LASSO | 0.4174 | 0.458 | 5 |
| Random Forest | 0 | 0.5477 | 19 |
for(i in 1:(nmodels-1)) {
dotchart(listr$CoefficientMat[i,],labels=colnames(listr$CoefficientMat),
xlab="Coefficients",main=paste("Coefficients - ",modnames[i]))
} # Line plot to the Predictions
par(xpd = TRUE, mar = c(2,2,2,7))
plot(1:n,listr$PredMat[,1],lty=1,type="l",
main="Predictions ordered by observation ID",xlab="Observation ID",ylab="Prediction",
ylim=c(min(listr$PredMat),max(listr$PredMat)), col=1)
for(i in 2:nmodels){
lines(1:n,listr$PredMat[,i],lty=i,col=i)
}
legend("topright",1,modnames,col=1:nmodels,lty=1:nmodels,cex=0.8,
inset =c(-0.18,0))
abline(v=92)Sur ce premier graphique, l’axe des x représente l’ordre dans lequel les mesures ont été effectuées. la droite verticale représente un changement dans le kit utilisé pour mesurer les expresisons des miRs. On peut constater que cela a un impact non négligeable sur les prédictions en cross-validation.
par(xpd = TRUE, mar = c(2,2,2,7))
idx <- order(classY)
plot(1:n,listr$PredMat[idx,1],lty=1,type="l",
main="Predictions ordered by class",xlab="Observation ID",ylab="Prediction",
ylim=c(min(listr$PredMat),max(listr$PredMat)), col=1, xaxt="n")
axis(side=1, at = 1:n, labels = idx, las= 2, cex.axis=0.8)
for(i in 2:nmodels){
lines(1:n,listr$PredMat[idx,i],lty=i,col=i)
}
legend("topright",1,modnames,col=1:nmodels,lty=1:nmodels,cex=0.8, inset =c(-0.18,0))
abline(v=57)L’axe des x représente un vecteur ordonné de la classe des observations. La ligne verticale délimite les témoins des malades.
# Scatter Plot matrix of predictions
pairs(listr$PredMat,main="Scatter Plot Matrix of predictions")# Courbes de ROC
cutoff <- seq(min(listr$PredMat)-0.1,max(listr$PredMat)+0.1,by=c(0.01))
nco <- length(cutoff)
truepositive <- matrix(nrow=nco,ncol=nmodels)
falsepositive <- matrix(nrow=nco,ncol=nmodels)
npos <- sum(classY==1)
nneg <- sum(classY==0)
listpos <- (classY==1)
listneg <- (classY==0)
for(i in 1:nco){
for(j in 1:nmodels) {
# Rate of true positives in the set of positive
truepositive[i,j] <- sum(listr$PredMat[listpos,j]>=cutoff[i])/npos
# Rate of false positives in the set of negatives
falsepositive[i,j] <- sum(listr$PredMat[listneg,j]>=cutoff[i])/nneg
}
}
plot(falsepositive[,1],truepositive[,1],type="l",main="ROC Curves Q-PCR Models",
xlab="False Positive Rate",ylab="True Positive Rate")
for(j in 2:nmodels){
lines(falsepositive[,j],truepositive[,j],type="l",col=j,lty=j)
}
legend(0.6,0.75,dimnames(listr$PredMat)[[2]],col=1:nmodels,lty=1:nmodels)# Recherche des seuils optimaux qui minimisent la distance au point (0,1)
cutoffopt=listr$PredMat[1,]
for(i in 1:nmodels){
distances2 <- falsepositive[,i]^2+(1-truepositive[,i])^2
cutoffopt[i] <- cutoff[which(distances2==min(distances2))][1]
}
pander(t(cutoffopt))| Logistique-SW | PLSR | O-PLS | Logistique-RIDGE | Logistique-LASSO |
|---|---|---|---|---|
| 0.5572 | 0.5472 | 0.6272 | 0.5472 | 0.5472 |
| Random Forest |
|---|
| 0.007205 |
# Dessin des cutoff sur un scatter plot
for(i in 1:nmodels){
plot(YC,listr$PredMat[,i],main=paste("Cut-offs - Method:",
dimnames(listr$PredMat)[[2]][i]),ylab="Prediction")
abline(h=cutoffopt[i])
} PredictClass <- matrix(nrow = n, ncol = nmodels)
dimnames(PredictClass) <- dimnames(listr$PredMat)
TrueClass <- rep(1,n)
for(i in 1:nmodels){
PredictClass[,i] <- listr$PredMat[,i]>=cutoffopt[i]
}
# Calcul les nombres de bien classés avec les seuils optimaux
ClassTable <- matrix("",nrow = nmodels * 3, ncol = 2)
dimnames(ClassTable)[[1]] <- rep(c("","Class0","Class1"),nmodels)
dimnames(ClassTable)[[2]] <- c("FALSE","TRUE")
dimnames(ClassTable)[[1]][((1:nmodels)*3)-2] <- modnames
for(i in 1:nmodels){
ClassTable[(i-1)*3+(2:3),] <- table(classY,PredictClass[,i])
}
cat("Classification Tables")## Classification Tables
pander(ClassTable)| FALSE | TRUE | |
|---|---|---|
| Logistique-SW | ||
| Class0 | 42 | 15 |
| Class1 | 12 | 51 |
| PLSR | ||
| Class0 | 39 | 18 |
| Class1 | 11 | 52 |
| O-PLS | ||
| Class0 | 38 | 19 |
| Class1 | 10 | 53 |
| Logistique-RIDGE | ||
| Class0 | 43 | 14 |
| Class1 | 13 | 50 |
| Logistique-LASSO | ||
| Class0 | 42 | 15 |
| Class1 | 16 | 47 |
| Random Forest | ||
| Class0 | 57 | 0 |
| Class1 | 0 | 63 |
# accuracy
pc.accuracy <- c()
for(i in 1:nmodels){
pc.accuracy[i] <- sum(diag(table(classY,PredictClass[,i])))/n
}
pc.accuracy <- as.matrix(pc.accuracy)
colnames(pc.accuracy) <- "pc.accuracy"
rownames(pc.accuracy) <- modnames
pander(pc.accuracy)| pc.accuracy | |
|---|---|
| Logistique-SW | 0.775 |
| PLSR | 0.7583 |
| O-PLS | 0.7583 |
| Logistique-RIDGE | 0.775 |
| Logistique-LASSO | 0.7417 |
| Random Forest | 1 |